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ABSTRACT 

A class of numerical dissipation models for central-difference schemes constructed with 
second- and fourth- difference terms is considered. The notion of matrix dissipation asso- 
ciated with upwind schemes is used to establish improved shock capturing capability for 
these models. In addition, conditions are given that guarantee that such dissipation models 
produce a TVD scheme. Appropriate switches for this type of model to ensure satisfaction 
of the TVD property are presented. Significant improvements in the accuracy of a central- 
difference scheme are demonstrated by computing both inviscid and viscous transonic airfoil 

flows. 
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I. Introduction 


Central difference type schemes are currently being applied on a regular basis m 
the solution of the Euler and Navier-Stokes equations. A numerical dissipation model is 
included in these schemes, and it plays a crucial role in the determination of their success. 
The form of the dissipation model is quite often a blending of second-difference and fourth- 
difference dissipation terms. The second-difference terms are used to prevent oscillations at 
shock waves, while the fourth-difference terms are important for stability and convergence 
to a steady state. There is a constant to be specified for each contribution. However, 
using the model in conjunction with appropriate numerical procedures, these constants 
can usually be selected and maintained for a fairly wide class of fluid dynamics problems. 
This dissipation model allows shock waves to be captured with smearing over three to four 

mesh cells. 

Even though these central-difference schemes have proven to be reasonably effective 
in many cases, there are still strong motivations for reducing the numerical dissipation 
being produced. For example, by appropriate reduction of the artificial dissipation, shock 
wave representation and boundary-layer definition (especially the wall shear stresses) can 
be improved on coarse meshes. Such improvements in accuracy are especially beneficial 
for complex three-dimensional flows, which can demand extensive computational effort. 
In addition, better estimates of the limit of infinitely fine mesh values of aerodynamic 
coefficients for flows with shocks can be obtained. Also, the standard model has difficulties 
in hypersonic flow. Finally, for some problems, the influence of numerical dissipation needs 
to be severely limited in certain smooth regions of a flow field (i.e., near the trailing edge 
of an airfoil), while still maintaining stability near discontinuities. This difficulty cannot 
generally be resolved by simply reducing the global constants in the dissipation model. 

One can appeal to ideas from upwind schemes to improve the dissipation model, 
especially in the vicinity of shock waves. Upwind algorithms utilize concepts from charac- 
teristic theory in order to determine the direction of spatial differencing. They have been 
extended to systems of conservation laws using such approaches as the flux-vector splitting 
of van Leer [l] and the approximate Riemann solver of Roe [2]. A fundamental feature 
of these schemes is the use of a matrix evaluation of the dissipation either in the implied 
or direct sense. In so doing, the dissipative terms of each discrete equation are scaled by 
the appropriate eigenvalues of the flux Jacobian matrices of the Euler equations, rather 
than the same eigenvalue as in the dissipation model employed with central-difference 
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schemes. Also, upwind schemes can be designed to have the total variation diminishing 
(TVD) property, which prevents the occurrence of spurious oscillations. The disadvantage 
of these schemes is that, in general, they increase the operational count for processing 
mesh points by about a factor of two over that required by central-difference schemes. 
One would certainly like to more closely imitate the highly desirable behavior of the up- 
wind algorithms near flow discontinuities, and at the same time, retain the more efficient 
central-difference scheme over significant portions of a flow field. In addition, one would 
like to have the high degree of numerical efficiency that has been achieved by combining 
a central-difference scheme with a Runge-Kutta time marching algorithm, which includes 
residual smoothing and multigrid acceleration techniques. 

The primary purpose of this paper is to construct a numerical dissipation model 
for a central-difference scheme that has both the properties of matrix dissipation and of 
TVD. As a starting point, we consider the elements of a widely used scalar dissipation 
model. Modifications of this model that facilitate accurate viscous flow computations 
are also examined. In the next section of the paper, the intimate connection between 
the formulation for an upwind scheme and a centered-difference scheme is presented, so 
as to establish a foundation for a matrix dissipation model. Then a theorem is proved 
that provides a simple sufficient condition to determine when a central-difference scheme 
with dissipation terms comprised of second and fourth differences is TVD. In the following 
section, appropriate flux limiter functions consistent with the central-difference dissipation 
model are discussed. A multistage time-stepping scheme used in applications is next briefly 
described. Finally, numerical results are shown to demonstrate the benefits of using the 
matrix dissipation model. Both inviscid and viscous transonic airfoil flows are computed. 

II. Scalar Dissipation Model 


The basic elements of the scalar dissipation model considered in this paper were first 
introduced by Jameson, Schmidt, and Turkel [3] in conjunction with Runge-Kutta explicit 
schemes. This dissipation model has been used by many investigators [4-6] to numerically 
solve the Euler equations for a wide range of fluid dynamic applications. The same type of 
dissipation model has been applied to alternating direction implicit (ADI) schemes [7] and 
LU factored implicit schemes [8]. Several modifications of the model have been investigated 
in [9] and [10] in order to improve it and make it suitable for obtaining accurate and efficient 
solutions of the Navier-Stokes equations. In this section, the basic model and important 
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modifications are briefly reviewed. 

Consider the Euler equations in the form 


W t + fx + 9y = °> 


( 2 . 1 ) 


where W is the four-component vector of conserved variables, and f,g are the flux vec- 
tors. The independent variables are time t and Cartesian coordinates (m). W 
transformed to arbitrary curvilinear coordinates i - £(*.») end 1 . 

obtain n\ 

{J-'W^ + Fz + G^ 0 , 

where J _1 is the inverse transformation Jacobian, and 

F = fVr, ~ 9 x rf> G = 9 x l~ fVt- 

In a cell-centered, finite-volume method, (2.2) is integrated over an elemental volume in 
the discretized computational domain, and J ‘ is identified an the volume of the cell. 
Equation (2.2) can also be written as 

J~ l W t + AWz + BWr, = 0, 

where A and B are the flux Jacobian matrices defined by A - dF/dW and B dG/ 

To advance the scheme in time we use a multistage scheme. A typical step of a 

Runge-Kutta approximation to (2.2) is 

W W =W (o)_ ak p. [o t F^ + D'.GV'-V _ AD \ , (2.3) 

where De and D v are spatial differencing operators, and AD represents the artificial 
dissipation terms. The dissipation terms are a blending of second and fourth differences. 

That is, . 

AD = (D\ + D*„-D\- Dj) W, ( 2 ' 4 > 


where 


D\W = V € [(A, + ,, y e^ ( .) A € ] W iJt 
D\W = V € A * V * A *] 


(2.5) 

( 2 . 6 ) 
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and A { , are the standard forward and backward difference operators respectively as- 
sociated with the { direction. The variable scaling factor A was originally chosen as 

*"+i,r “ 2 + + (*e).,, + Mi+i,;] • (2.7) 

where A f and A„ are proportional to the largest eigenvalues of the matrices A and B. The 
scaled spectral radii A^ and A^ are given by 

A € = l^-^l + c^ + ya, 

= \vxt - uy f | + cyjx\ + y|, 

where « and V are the Cartesian velocity components, and c is the speed of sound The 
coefficients *(*> and «<«> are adapted to the flow and are defined as follows: 

~ ^ ^ ^,+i 1 y,i/, + 2 i y), (2.8) 


‘ ,y \*+ij + 2Pij + Pi-ij\' (2-9) 

e ‘+l.f = maX [°> - &’*,,•)] ' (2.10) 

where p is the pressure, and the quantities *<»> and ,=M are constants to be specified The 
operators in (2.4) for the , direction am defined in a similar manner. 

The second-difference dissipation term is nonlinear. Its purpose is to introduce an 
entropy-1, ke condition and to suppress oscillations in the neighborhood of shocks. This 
term is small in the smooth portion of the flow field. The fourth-difference dissipation term 
is basically linear and is included to damp high-frequency modes and allow the scheme to 

approach a steady state. Only this term affects the linear stability of the scheme. Near 
shocks it is reduced to zero. 

The isotropic scaling factor of the original dissipation model as given in (2 7) is gen- 
erally satisfactory for inviscid flow problems when typical inviscid flow meshes (i e cell 
aspect ratio 0(1)) are used. The factor can produce too much numercial dissipation in the 
cases of meshes with high aspect ratio cells. This is also an important consideration for 
high Reynolds number viscous flows, where a mesh providing appropriate spatial resolution 
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can have cell aspect ratios O(10 3 ). In [9] and [10] this difficulty is remedied by replacing 
the factor of (2.7) with the anisotropic one 

= 2 [fidiJ + (^)t+lj] 

where 

= M r ) ( A f 

= 1 + 'ij 0 < f < 1, 

and r = A n /A^. In the normal direction, one defines 

Alternatives to the switching function presented in (2.9) have been investigated. Cau- 
tion must be exercised in the selection of a switching variable. If a quantity with the same 
functional dependence as entropy (i.e., p/p' 1 ) is used, sharper shocks can be captured in 
viscous transonic flows. However, such a choice can result in a loss of accuracy for the 
surface shear stress, due to the significant variation in the entropy-type variable across the 
boundary layer. This difficulty can be removed by simply multiplying the scaling factor 
by a function of the local Mach number of the flow. An acceptable modifying function has 
proven to be (Ml/Mqo)^, for some (3 > 1, where M L is the local Mach number, and M re f 
is a reference Mach number (i.e., free-stream value for external flows). It is important 
to note that the entropy-type function is generally not satisfactory for inviscid flows. In 
addition, one can consider the sum of two switches, one depending on pressure and the 
other on temperature so that all thermodynamic changes are taken into account. When 
introducing the matrix-valued dissipation it will be possible to use separate switches for 
different characteristic variables. 

The treatment of the artificial dissipation must be modified at the boundaries of the 
physical domain. In the case of the fourth-difference dissipation the standard five point 
difference stencil must be replaced at the first two interior mesh cells. This means that 
one-sided or one-sided biased stencils are used at these cells. The dissipative character 
of the artificial terms is important because it influences both stability and accuracy. For 
example, if the dissipation is too large at a solid boundary, an artificial boundary layer is 
created in an inviscid flow, and the effective Reynolds number for a viscous flow is altered. 
To improve accuracy at the wall boundaries of viscous flows, where gradients are steep 
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due to physical boundary layers, the usual fourth-difference stencils are changed in this 
dissipation model. 

Let the total dissipation for a mesh cell, in the direction represented by the index j , 
be denoted by dj. For simplicity assume that Ae( 4 ) = 1. Then, 

dj = d/y + 1 — dfj_ l. 

where the dissipative flux 

df j+i = (AW) 1+i - 2(AW) j+i + 


and thus 


d > = - 3(A W) j+i + 3(AW) j _ i - (AW); 

-f (2-11) 

with the index * for A W suppressed for convenience. Consider the first two interior cells 
adjacent to a solid boundary, as depicted in Figure 1. If 

(AW) l. = (AW) a = (AW)s 

(2.12) 

then (2.11) gives 


d 2 = W 4 - 2 W 3 + W 2 , 

(2.13) 

d 3 = W 5 - 4 W 4 + 5W 3 - 2W 2 . 

(2.14) 

These boundary stencils are fairly standard ones, and they result in a nonpositive definite 
dissipation matrix for the system of difference equations [7]. An alternative form, which has 
reduced the sensitivity to solid surface normal mesh spacing for turbulent flow calculations 
without compromising stability or convergence, is given by 

(AW)i = 2(AW)| — (AW)| 

(2.15) 

and 


d 2 = W 4 - 3W 3 + 3W 2 - W u 

(2.16) 

d 3 = W 5 - 4 W 4 + 6W 3 - 4W 2 + Wi. 

(2.17) 


This boundary condition is advantageous if the mesh is fine enough to adequately represent 
the laminar sublayer region of the boundary layer (i.e., at least two mesh points are inside 
the sublayer). For coarse meshes this treatment can be less accurate than the zeroth-order 
extrapolation of (2.12). 
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III. The Upwind Connection 

Upwind schemes for solving hyperbolic systems of conservation laws (i.e., Euler equa- 
tions of gas dynamics) rely upon characteristic theory to determine the direction of prop- 
agation of information, and thus, the direction required for one-sided differencing approxi- 
mations of the spatial derivatives. With such schemes shock waves can be capture 
oscillations. Thus, a successful artificial dissipation model for a central-difference scheme 
should imitate an upwind scheme in the neighborhood of shocks. We now r 

connection between these two types of schemes. 

Consider the one-dimensional scalar wave equation 

Ut + au x — 0 

with a constant. The first-order upwind scheme can be written as 

“<° ( 3 . 1 ) 

U n + 1 =u j -a—- < V ' 

1 Al Ui-uy-„ o>0, 

where all discrete quantities are evaluated at time level »At unless otherwise denoted. The 
scheme of (3.1) can be rewritten as 

(3.2) 


“” +1 = - “2S (Uj+1 _ “*- l) + |o| 5S (Uj+1 ” 2u ’ + 

which now contains a central-difference term and a second-difference dissipation term. Now 

consider the system . 

u t + Au x = 0, W 

where » is a N-component vector. The system case can be converted to a scalar one by 
diagonalizing the N x N matrix A with a similarity transformation 

A = AT, 

where the columns of T are the right eigenvectors of A. After diagonalizing (3.3) and 
applying the scheme of (3.2), the first-order upwind scheme is given by 

“ y +1 = «» ' A 2^ ( “ i+1 " + |A| ^ (U ’ +1 “ 2u ’ + u ’~ l) ' <3 ' 4) 


where 


A = Diag [|Ai| • • • |Ajv|] ■ 


\A\ = T|A|T _1 , 
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The generalization to a system of conservation laws is straightforward; namely, 


«* + /* = 0 


with / being a N-component flux vector, and 




^ ^ + 2Ax [K+t ( u i+i ~ «/) “ (uy — tiy-i)] , (3.5) 


where the Jacobian matrix A = df/du, and \A\ is defined as for (3.4). The matrix 
| Ay+ *l can be computed as an arithmetic average or a Roe average. For transonic steady 
flows the differences are negligible; therefore, we use the simpler arithmetic average. For 
hypersonic flows Yee [11] found that the Roe average yields better results. For time- 
dependent problems the Roe average also seems to give slightly better results. 


IV. Matrix Dissipation Model 


We now extend the scheme given in (3.5) to the two-dimensional equations of fluid 
dynamics. In particular, consider the transformed Euler equations of (2.2) with the Runge- 
Kutta scheme of (2.3). The necessary modification to the contributions for the £ direction 
of the artificial dissipation term defined by (2.4) is to substitute \A\ for the eigenvalue 
scaling factor, A, in (2.5) and (2.6). For the r, direction, £ and \A\ are replaced by r, and 
\ B \, respectively. We next define explicitly the form for the matrix \A\. Let 


with 


Then, 


where 


A Diag [Ax A2 A3 A3] 


A i=q+<Ja\+alc, \ 2 = q — a\ + a.\ c, A 3 = «, 

a i = o 2 = £y, q = aju + a^v. 

w - !*•!'+ - |As|) fcb, + ^E. 

' 1 / L c 2 a \ + a\ 


+ 


l^ll ~ |^2 1 


V T?+a|c ) |g3 + h - 1 > g< l- 


E , = 


<t> 

u<f> 

v<f> 

H<j> 


—u 

-u 2 

—uv 

-uH 


1 1 

u 

V 2 V 

vH H 


—v 

—uv 

.2 
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£'2 — 


E$ — 


E 4 = 


■ 0 

0 

0 

o- 


-a iq 

a\ 

a io 2 

0 


—a?q 

d 2 d\ 

a 2 

0 

9 

L-g 2 

qa x 

qd 2 

0. 


' -q 

d\ 

0*2 

(T 


—uq 

ua i 

Ud 2 

0 


—vq 

Vdi 

Vd2 

0 


.-Hq 

Hdi 

Hd 2 

0. 


■ 0 

0 

0 

0 

1 

ai<f> 

-ai u 

—ait; 

a x 

d2<t> 

— d2U 

— d 2 V 

o 2 

- q(f> 

—qu 

-qv 
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H is the total enthalpy, and <p = - tv — * t 

row of E, is a scalar times the first row) for any X„ X,. and X, an arbitrary vector, « 

be multiplied by |A| very quickly. That is, we calculate |a ;H | («f+i «f) ( 

(1 2] ) rather than calculate |a #+ j| and multiply a matrix times a vector. The matnx |B| 

is computed in the same way as \A\ by simply replacing £ with r). 

* t, ol i i,as riven above. Near stagnation points A 3 
In practice one cannot choose A x , A 2 ,A 3 as given a 

p . .. a \ anoroach zero. A zero artificial viscosity 

approaches zero while near sonic lines A x or A 2 approach zero. 

would create numerical difficulties. Hence, we limit these values as 

|\,| = max(|X 1 |,V' n p(A)), p(A) = 111 + c\/«? + “I 

(A, | = max(|X 2 |,V n p(A)), |*»| = max(|X 3 |, V t p(A% 

where the linear eigenvalue X 3 can be limited differently than the nonlinear eigenvalues 
The parameters and Vr have been determined numerical, y. Vanous .values b«b« 
evaluated by comparing their corresponding computed solutions on the basts of 
ing . 1) Sharpness of shock waves captured (without producing oscillations), 2) Convergence 
rate of numerical scheme. A good choice for V n and Vr is between 0.2 and 0 3. 

We have thus far replaced X j+ ., in (2.5) and (2.6) by a matnx while leaving ,h 
limiters e< 2 > and d 4 ' as scalars. One can also introduce e 2 and e into e tagon 
matrix A. This allows different limiters to be chosen for different characteristic variably 
For example, the limiter may be based on pressure for the nonlinear waves. However, the 
pressure is smooth through a contact discontinuity. Hence, a switch based on temperature 
may be more appropriate for the linear wave. One could also use different mesh scaling , 
m for the linear and nonlinear waves. Also, a smoother cutoff (13] can be mtroduc . 
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V. The TVD Property 

Consider the one-dimensional scalar conservation law 

+ = (5.1) 

where 

-oo < x < oo, t > 0. 

Let v(t) = be the approximate solution of (5.1) and consider the semidiscrete 

equation 

di Vj W + 2A^ /j+1 ~ = 2Ki Av 3+k - 

*<«> r D i 

Ax v 3+h (5.2) 

with 

At V+* = ( Av W = ®i+i W - * v(<); 

A 3 is a third-difference operator defined as 


A ( A V }j+?~ V i+ 2(0 3wy +1 (t) 4- 3»y(t) — Uj-jJt). 

The terms on the right hand aide of (5.2) represent second- and fourth-difference numerical 
dissipation terms, with /c< 4 ) a constant. Define 


s >+* = s « n . 

where sgn represents the signum function. Then, a modification of the theorem of Tadmor 
[14] is given by the following. 

Theorem: The semidiscrete scheme of (5.2) is TVD if 


(a) 


and 


[’ 


+ 5 y+f)] Qj+l > «,■+*(. 


'i+W+i ~ s 3-k) 


a/ ; 




*' Au y+* 


0 >) 
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Proof: Shifting the indices hy one in (5.2) and subtracting (5.2) from the resu.ting equation, 


we obtain 


+ sj [ A/,+ * Afi '^ , 

= jij [0,- + *Av, + . - 2Q i+i A « i+i +Qj-l A«i-i] 

_ ^ [ E) . + , A 3 o )+ j - 2ff j+ . A 3 » J+ . + 


(5.3) 


Multiply (5.3) by s j+ . and sum over 


all j. Note that s y+ i. = ±1, so s 2 j+ i - 1, and 


;j _i Av 1+ i = Av /+£ 


S j+ L *V j+ 


We then obtain 


— V |Av 
dt 1 


A/j+i 1 AUj+^l 


y+ 5 1 


Av yH 


2Ai 


Av y+^ 


We stress that the last term in 
frequencies and accelerate convergence to 


5 >,+* (s )+ i - 25 , +i + »>-*) «>+ 

i 

3 

(5.4) will not help for TVD. Its purpose is to eliminate high 

' , J 1 * . _ X rkn 


(5.4) 


a steady state. Hence, we want this contribution 


to 


be zero. This can be accomplished if we demand either 


s j+| “ 2s i+k +s i-± 


= 0 


or 


*,h = 0 ’ 


(5.5) 


(5.6) 


i.e., condition (b). We are left with 


j i v' r t 

J t ( TV ) = ~2Ax ^ [- s y+*( 5 y+§ 8 i-tf Av j+h 

- s j+ ^{s j+ 1 - 2s y+ i. + s ,-^) < 5y+5 


Av y+* 


(5.7) 


where TV denotes the total variation as given by 



tie? " 8 “. b T ^ COnditi °” ‘ hat the tOU ' not -crease is that 

term of (5.7) in brackets must be positive. This means that 


s >+ir ( 5 y+| 2s y+* + *y-*) Qj+^ > s j+k (« y+| - a ._ h ) ^±£. ( 5 . 8 ) 

y+i 

Since a’ +J = 1, (5.8) is equivalent to condition (a) of the theorem. Defining 

Xy = 1 -*y-j«y + j 

(5.8) can be written in the form given by Tadmor; namely, 


(xy + Xy+i )Q i+ t > (xj - Xj+i) ~ 3+ ^ 


Aw, + 1 

j+i 


Several remarks concerning this theorem are in order. 

Remarks 

(i) When no extrema are present locally \ ... 

.1 .1 . ' J 5 y+f)’ both conditions of 

he theorem are trivially satisfied for all Q , , and R , pw c ,v u 
n — r\( \ \ c nd K i+±- For suc h regions we want 

i+i ~ (Ax) for second-order accuracy and R j+ , to be chosen so that high fre- 
quencies are damped. 

(“) If Sj ~* = *’+! = ~*>+h ( a local oscillation at x j+ x ), we require 

<?>+£> 0, R jH= 0. 

(" ") y 5 y+i s y+i- ( a loc&l extremum 


■ ? y+ 2 

Qj+h- > 


at Xy +1 ), 


I^y+l 

* “ A ”y + ^ ’ 


R-i+L =0. 


(iv) If«y_^ 5 j+» — s J + 3 (a local extremum at xy), 


Qi+i > 


' j+ * ~ Aw,, , * 


R"i4~ i* = 0. 


It follo»s from these inequalities that Q J+ , can be negative. As far as total variation is 
concerned, central d, (Terences are not nondissipative. That is, they can either increase or 
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decrease the total variation. In cases where centra, differences decrease the total varia- 
tion, Q, + . can be negative. For systems of equations, especially m multuI.mens.ons, .s 
behavior can sometimes lead to difficulties. Hence, we demand the stronger cond.fon a 


<3f+ i a 

For systems this condition is replaced by 




Av 


’}+? 


Qj+k - K+*|’ 

where A = df/du, and the average at i + 1 is constructed by the technique of Roe [21. 
VI. Flux Limiters 

In this section switching functions are introduced that force the scheme to automat- 
ically satisfy the inequalities presented in Section V. These switches are required to be 
smooth so that limit cycles are not experienced when marching in time to obtam a steady- 
state solution. The discrete switching functions are defined as 


Q j+k = tl>-\aj+L | with 0 < 0 < 1, 


(6.1) 




a,-, i 


where = 1 near extrema, so that the inequalities are satisfied, and 4> = O(Ax) in smooth 
regions of the How field. Conversely, T = 0 at extrema, and T = 1 in the smooth regions. 
The functions and I' of (6.1) can be defined in terms of a limiter function. Let 


r = 


_ Vj - Vj - 1 _ ^ v j- k 


Vj + 1 


— Vv 


Av.j » 


The van Leer flux limiter is given by 


= iTFi “ \ 


v j+k 


-rr-, r > 0 

1+r ’ 


^ 0, r < 0. 

From Sweby [14] it is straightforward to see that for any flux limiter 

V>j( r ) = 1 “ 
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Since we want ^y(r) to be positive, this relation is redefined as 

^;'( r ) = |1 “ <P}(r)\. 

Then, for the van Leer limiter, 


1 - r 

i 

i + M 

i 


V-y(r) = 

so 0 < ipj < 1 and define 

't’i+k = max (V’y, V’y+i)- 

We now show that the inequality (5.8) is satisfied with this By the first remark of 
Section V, the inequality is satisfied when = s,. + , = a,. + } , and so in this case we 

only need Q — O(Ax). Since r > 0 in this 

V-y(r) = 


case, 

1 — r| 


1 + r 


In addition, for smooth regions of the flow field, 

r = 1 + O(Ax), 


and thus 


= °(Ax), Qy + i = O(As). 

Next, consider the cases,.. = -s, + , . This implies r, < 0, and so, ^ 


= 1. Moreover, 


5 y-^ 


' s y+i => tfy = 1 =J> 0y + ! = 1. 


Similarly, if s J + 3 = 


— s 


y+j> 


s y+^ ~ 5 y+| V’y+i = 1 => = l. 

It also follows, using (6.1), that the inequalities in (ii) and (iii) of the remarks of Section V 

are satafied. Also, for these two cases, setting T, + j = 1 -* /+ , guarantees that r . . . = 0 
So ■ 7+5 ' 

r y+§- = 1 -max(^y,0y +1 ). 
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The ffux limiters * and F can be connected to those used in upwind schemes. Sweby [15] 
considers an upwind Lax-Wendroff scheme. In particular, for the one-dimensiona wave 

equation 

Ut ~h Q>u x = 0, 

the numerical solution is obtained with 

I 1/(1 - v) 


where 


u" +1 = Uj - v{uj - «y- 1) - A - 


i/ = a^, 0 < <p{r) < 2, 

Ax 


-<Pj(Uj+l - Uj) 




(6.2) 


and the backward difference operator A_ is defined by 


A -Uj = Uj - tty— i. 


If (6.2) is rewritten as a central-difference scheme, then 

u n +1 = _ t {uj+l - + £[(1 - W) (•,« - -i) - d - W-‘) - U ’- ,)1 

3 2 


+ - 2 [v>y(tty+i - tty) - <Py-i( u y _ “j- 1 )] • 
2 


(6.3) 


By dropping the term in (6.3) and changing to a semidiscrete formulation, we get 

+ — [(1 - < pM «,+. - « i ) - (1 - *>»-»)(“>• - “ t - 11 • 

2Af 


If 

<Pj = 1 “ V»y+i> 

then the second-difference dissipation term has the same form as the one presented in 

Section V. ^ , • j 

To complete the connection between the limiters for the central-difference and upwmd 

schemes, we compare their behavior. For the central-difference case, 

o<V>< i, o< ^ 

whereas for the upwind case, 0 < < 2. Furthermore, for the central-difference limiter, 




(1) = *(r), if (i) = «’(’■) 
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means that it does not matter if Au, + j > Alty.j or Au, + j. < (i.e., the sign of a 

oes not affect the scheme), as opposed to the upwind limiter where 


<p( 


'(r) fl\ 


In Figure 2 the Sweby diagram for the upwind van Leer flux limiter and the central- 
difference van Leer flux limiter is shown. For r < 1, the limiters are the same. For r > 1 
the upwind version continues until p = 2, while the central-difference version returns to 
zero. In Sweby s paper the flux limiters are not allowed to decrease when r > 1. However 
there is no difficulty from such behavior. In both cases p(r) = 0 when r < 0 so that the’ 
limiter is turned off for extrema. We note that other limiters besides that of van Leer can 
be used to get switches for central-difference schemes. 

For systems of equations we use a scalar limiter. Using the matrix form of the dissi- 
pation it is easy to implement different limiters for different characteristic variables. 

In terms of the pressure the switch becomes 


V»y = 


1 - r 
1 + 1 rl’ 


r = 


Ap y _ 


or 


with 


= jpy+i -Pj\- ( Py - Py-i) Sgn (p j+l - P-) 
IPj+i ~Pj | + |py — Pj-i \ + e 


a > 0 


a > 0 


ti+k = max (^y, i>j+ 1), 

and c = OfAx 2 ) to prevent a zero denominator for constant pressure regions. Consider 
also the wave speed a < 0. Then, 


tjj- = Pj+i - 2 py + pj_i 

\Pj+ 1 ~ Pj\ + by - Py-i | + e 

We use a conservative approach and take 


sgn(p y+ i - Pj ) t a > 0 
s g n (Py-i - Py), a< 0. 


xb- = IPj+i ~ 2py + Py-i| 

J |Py+i -Pyl + |Py-py-i| + e' ( 6 - 4 ) 

Notice that this switch is very similar to (2.9) for the original dissipation model of the 
Runge-Kutta scheme. There is only a minor change in the denominator. However, with 
this change and the factor 1/2 in front of the second-difference dissipation term, the scalar 
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equation becomes amt-order upwind near shock. In the case of the original * we find 
that « .05 near shock waves in transonic Bows. One may require different paramet 

for the Runge-Kutta scheme to ensure stability. 

We now no ionger have a free parameter for the second-difference a, 

usually use T = 1 - H so that the fourth-difference dissipation is cut off for ^ > 1/2. 
only free parameter is the coefficient *<*> of the fourth-difference term. 

VII. Numerical Algorithm 

The majority of the numerical results presented in this paper were obtained with a 
Navier-Stokes code developed by the authors, which is based on the explicit mult, stage 
time-stepping schemes of [3] and [16]. This class of schemes » current y in w. esp 
use for solving the Euler equations. In references |17-20], these schemes were extende 
to allow solution of the compressible Navier-Stokes equations. Significant 
in numerical efficiency were introduced in [9], [10], and [21], In the code o Swanson and 
Turkel a cell-centered, finite-volume method is employed to obtam “ntere ype . 
approximations for the flow equations. Such a method provides flexrb.hty m trea g 
arbitrary geometries and different grid topologies, since no special treatment is required m 
the vicinity of singular points or lines. The scheme is second-order accurate in space or 
sufficiently smooth meshes (see [21-22] for definition of sufficiently smooth). A modiBe 
five stage Runge-Kutta scheme is used for the time integration to obtam a steady- 
solution. There is a weighted evaluation of the artificial dissipation terms on the first 
third, and fifth stages, which provides a large parabolic stability limit. e P Y 
viscous terms are computed on the first stage only and frozen for the remaining ones; is 
does not compromise the stability characteristics of the scheme. The spatial and temporal 
differencing are decoupled. Thus, the numerical algorithm is independent o time s ep 
and amenable to steady-state convergence acceleration techniques. These methods inclu e 
local time stepping (a preconditioning for the system of difference equations), variable 
coefficient implicit residual smoothing, and multigrid. Implicit residual smoothing is jus 
a mathematical step, applied at each stage of the explicit scheme, to extend the local 
stability range. The multigrid method involves cycling through a sequence of successive y 
coarser grids and relying upon effective high frequency damping for rapid removal of 

errors in the fine grid solution. 
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VIII. Applications 


Two airfoil flow problems are considered here to demonstrate the benefits of using a 
central-difference scheme with a matrix numerical dissipation model. The first problem 
concerns inviscid flow over an NACA 0012 airfoil. At a free-stream Mach number (M„) of 
0.8 and angle of attack (a) of 1.25 degrees, a fairly strong transonic shock wave occurs on 
e upper surface, while a much weaker one appears on the lower surface. On representative 
inviscid meshes, schemes based on central differencing capture the upper surface shock 
reasonably well, but they smear significantly the lower surface shock. The second problem 
involves transonic turbulent flow over a RAE 2822 airfoil. In this case the free-stream 
Mach number is 0.73, the angle of attack is 2.79 degrees, and the Reynolds number based 
on c ord (Reoo) is 6.5 x 10*. For such transonic viscous flows small differences in the shock 
strength can result in noticeable changes in the lift and drag of the airfoil. Thus, both 

o t ese problems can provide a reasonable measure of the performance of the artificial 
dissipation model, especially near shock waves. 

A C-type mesh consisting of 224 cells around the airfoil (160 cells on the airfoil) 
and 32 cells normal to the airfoil was used for the first problem. The outer boundary 
of the finite domain was placed 20 chords away from the airfoil, so as to not produce 
significant effects on the solutions (see (23|). The normal mesh spacing at the surface was 
approximately 0.01 chords. The mesh was clustered at the leading and trailing edges of the 
airfoil m order to improve resolution of the steep gradients occurring in these regions In 
Figure 3 the surface pressure distributions computed with the scalar and matrix dissipation 
mo es are shown. There is a discernible improvement in the sharpness of the shock 
waves using the matrix model, especiafiy for the one on the lower airfoil surface. The 
lift and drag coefficients obtained using the two models, along with the values from a 
high density mesh calculation, are compared in Table I. The values associated with the 
matrix model nearly match those corresponding to the fine mesh solution. The convergence 
ist ones m terms of number of multigrid cycles am also displayed in Figure 3. Convergence 

“ meaSUred by the loBarithm ° f thc «** mean square of the residual of the continuity 
equation. The convergence rates obtained with the scalar and matrix models are 0 819 

and 0.888, respectively. The matrix model result required about 55 cpu sec on a Cray II 
computer for 125 multigrid cycles with the solution grid This processing time represents 
s percent increase in the time needed compared with the scalar dissipation model. It 
may be possible to improve the convergence rate with the matrix model by altering the 
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. a //\t- fViA tvne of implicit residual smoothing 

coefficients of the time-stepping scheme and/ V 

operator employed. 

The following set of C-type meshes were used in solving ‘he -cond problem with 

• a i • ifiO x 32 with 128 cells on the airfoil, 2) 32 

thp two dissipation models: 1) loU x az wmi t 

256 cells on the airfoil, 3) 640 X 128 with 512 cells on the airfoil. Each successive y coars 
mesh was generated by eliminating every other mesh line in both coordinate irec ions 
r air mesh. Again the outer hou„aary of the domain was .ocated . .ds from t 
. . -l The normal spacing at the surface for the Bnest mesh was 7.5 x 10 chords. 

4 a comparison is made between the experimental data of | M ] and the predicted 
distributions of pressure, upper surface skin friction, and upper surface boundary ^ 
displacement thickness obtained on the 160 x 32 mesh. The skin fnct.on -ment.the 
wail shear stress nondimensionalized by the dynamic pressure at the edg 
"ocity proa.es (scaied hy the boundary-, ayer edge velocity «.) at two stations on 
ripe, a JoU surface are a, so presented in Figure 4. There is marked improvement n 
the results obtained with the matrix mode,, as verified hy their —b e agreemen with 
the data. The convergence histories given in Figure 4 indicate that the final residua, lev 
‘"somewhat higher with the matrix mode,. However, the termina, rate of convergence is 
nearly the same with both dissipation models. Finally, Figures 5 and 6 show the pressure 
and skin-friction distributions computed on the coarse mesh with second-order upwin an 
“third-order” upwind biased forms of Roe's scheme |2], which has been shown to have low 
levels of dissipation. The pressures obtained with the “third-order” form are « 

those calculated with the matrix model. The upwind skin-friction solutions exhibit slightly 

better agreement with the experimental data upstream of the shock. 

Numerical results for the two finer meshes are presented in Figures 7 and 8. For the 
320 X 64 mesh, a slightly stronger shock is predicted using the matrix issipa ion m 
Otherwise, the solutions determined with the different dissipation models * 
same The character of the convergence behavior with the two models is about 
as for the coarse mesh. On the finest mesh, the application of the models gives essentially 
the same results. In Figure 9 the variation of the computed lift and tota rag coe jn 
with the reciprocal of the number of mesh cells is plotted for each dissipation m . 

The curves corresponding to the matrix model are nearly linear. A linear curve indicates 
second-order accuracy for the full range of meshes being considered. The numencal va 
for the components of the drag coefficient (form and friction contributions) along with 
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hft ccfflcent are given for all the viscous computation* in Table II. Notice that the value* 
obtamed with the matrix model on the 320 x 64 grid are very close to those calculated 

‘ hC SCallr m ° deI °" thc 640 * 128 S'*- Such a reduction in the mesh required for 
acceptable accuracy is highly desirable, especially for thre«Iimensional Bow problems 

ome apphcatmns to three-dimensional viscous flows are presented in Turkel and Vatsa 

! : ' Sh ° Uld be emphaSized that ““ — prediction of lift and drag is very important 

e esign of aircraft, and thus, a good estimate of these quantities for an infinitely fine 
mesh is needed. 

In [26| it is shown that the TVD switch (6.4) allows converged solutions for hypersonic 
flow where the standard switch (2.9) does not. 

IX. Concluding Remarks 

We have thus shown that the second-difference artificial dissipation is equivalent to 
using a flux I, miter, and hence, central-difference schemes are not any more “artificial” 
than upwind schemes. The central-difference scheme is slightly more dissipative for two 
reasons. First, - max(i/i ; , 0 J+] ) while for an upwind scheme ^ +1 is equal to either 

,! OI ^ +1 ’ depen<Ung ° n the dire<:ti °" <>f *!>« win* Second, we insist that be positive 
(■•e., (Py < 1) while upwind limiters allow > 1 (i.e., in some cases we can have negative 
viscosity but still be TVD). However, to compensate for this slight increase in dissipation 
central-difference schemes are simpler to program and require less computer time per time 
step, and work well with multigrid acceleration techniques. 

th f In thc centraMi ' rerenc ' schcmcs w ^ ««*««««, ^ 

the fourth-difference dissipation. This dissipation is needed to approach a steady state 

and has nothing to do with TVD properties. In fact the fourth-difference contribution is 

set equal to aero near local extrema. For time-dependent flows one can set this dissipation 

y Zer °- ° n the ° ther hMd TVD P ro Perties do not necessarily imply a rapid 
convergence to the steady state. P 

In summary, a formulation for a numerical dissipation model that makes a central- 

' T“ ClOSely reSemWe “ u <> wM "uur flow discontinuities has been 

described. A theorem has been proven that gives sufficient requirements for this type of 

dissipation mode, to satisfy the TVD property for a scalar equation. Flux limiter functions 
ave been presented for this form of dissipation model. For a system of equations a matrix- 
valued dissipation is introduced. Solutions of the Euler and Navier-Stokes equations for 
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airfoil Bows have been obtained using .he matrix dissipation mode,. These results have 
demonstrated noticeabie improvements in accuracy in smooth regions . the flow fle d 
as well as near shock waves. There is a 15 percent increase in computational 
explicit multistage schemes when this model is used. However, the improved accuracy has 
permitted a signiflcan, reduction in the number of mesh points reared, uch behavior* 
I scheme can have a dramatic effect on the necessary mesh size for three-dimensional flow 

CalC “,' it is important to emphasize the different principal objectives associated with 
matrix-valued dissipation and the TVD switch. The purpose of the matrix form o e 
numerical dissipation mode, is to apply the appropriate scaling of the diss.pat.on in each 
flow equation, yielding a reduction in the amount of dissipation being introduce an 
improved accuracy. The TVD switch plays a somewhat opposite role in that it causes 
more dissipation to be added in order to prevent overshoots, and thus, allows convergence 
and provides robustness in solving high speed flow problems. The combination of the two 
should give more dissipation near shocks and less dissipation in smooth regions, hence 

giving better accuracy in all regions. 
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Table I 

Lift and drag coefficients for NACA 0012 airfoil. 
Moo = 0.8, a = 1.25° 


Case 

ci 

Cd 

Scalar dissipation model 
(224 x 32 mesh) 

0.3628 

0.0231 

Matrix dissipation model 
(224 x 32 mesh) 

0.3591 

0.0227 

Scalar dissipation model 
(640 x 64 mesh) 

0.3577 

0.0228 


ci - lift coefficient 
c d - drag coefficient 
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Table II 

Lift and drag coefficients for RAE 2822 airfoil, 
Moo = 0.73, a = 2.79”, Re„ = 6.5 x 10 6 


/IJccirtat.inn model 



160 x 32 
320 x 64 
640 x 128 

0.8296 

0.8514 

0.8588 

0.0124 

0.0123 

0.0124 

0.0050 

0.0055 

0.0055 

0.0174 

0.0178 

0.0179 

Second-order upwind 

160 x 32 

0.8176 

0.0119 

0.0054 

0.0173 

“Third-order” upwind biased 

160 x 32 

0.8220 

0.0124 

0.0054 

0.0178 


ci - lift coefficient 
Cd - pressure drag coefficient 
c df - friction drag coefficient 
C j - total drag coefficient 
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Figure 2 Sweby diagram 
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res,lIts for inviscid flow over NACA 0012 airfoil 
(a) Surface pressure coefficient, (b) Convergence history 
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Figure 4 Central-difference scheme results for turbulent 

(160 X 32 mesh, Moo = 0.73, a = 2.79 degrees, Ra» . - x 10 ). 

(a) Surface pressure coefficient, (b) Convergence history 
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( a ) ( b ) 


fl'eox 32 m X°Mt°= artl = h ™9 dTgrli R^ter: fo^-° Ver ““ 2822 airfoil 
(a) Surface pressure coefficient, (b) Upper surface skin-friction coefficient 
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x/c 



( a ) 


(b) 


Figure 6 “Third-order” Roe scheme results for turbulent flow over RAE 2822 airfoil 

(160 x 32 mesh, M«, = 0.73, a = 2.79 degrees - «; 5 ^ )• 

(a) Surface pressure coefficient, (b) Upper surface skin-fnction coefficien 
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( a ) ( b ) 


/of n Rd , t 1 * ™* scheme results for turbulent flow over RAE 2822 airfoil 
( 4 mesh, Moo — 0.73, a = 2.79 degrees, Re^ = 6.5 x 10 6 )- 

(aj Surface pressure coefficient, (b) Convergence history 
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□ EXPERIMENT, COOK 



Figure 7c - 7d Upper surface skin-friction coefficient (c,) and boundary-layer displacement 
thickness (320 X 64 mesh) 
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Figure 8 Central-difference scheme results tor turbulent flow over RAE 2822 airfoil 
(640 x 128 mesh, = 0.73, a = 2.79 degrees, Re,*, = 6.5 x 10 ). 

(a) Surface pressure coefficient, (b) Convergence history 


37 


Figure 8c 
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Upper surface skin-friction coefficient (640 x 128 mesh) 
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Figure 9b Variation of total drag coefficl 
(RAE 2822 airfoil, = 0.73, a = 2.79 de 
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